Tarea 2: Frequentist Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega: 08/10/2024

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Enlaces a videos de las clases:

Documentación:

Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Para realizar plots
library(ggplot2)
library(plotly)
## 
## Adjuntando el paquete: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## 
## Adjuntando el paquete: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

Elimine eval=FALSE para ejecutar las celdas.

Pregunta 1: Estimadores.

# Función para generar los estimadores
calcular_estimador <- function(n, sigma) {
  # Simulación de las X_i ~ Bernoulli(0.5)
  X <- rbinom(n, 1, 0.5)
  # Ruido epsilon_sigma ~ N(0, sigma)
  epsilon_sigma <- rnorm(1, mean = 0, sd = sigma)
  # Estimador
  p_gorro <- epsilon_sigma + (1/n) * sum(X)
  return(p_gorro)
}
# sigma = 0
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_0 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())
for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 0)
  
  # Almacenar los resultados
  resultados_0 <- rbind(resultados_0, data.frame(n = n, sigma = 0, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_0, aes(x = n, y = p_gorro)) +
  geom_line(color = "blue") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 0)",
       x = "n",
       y = expression(hat(p)[n, 0])) +
  theme_minimal()

# Sesgo para sigma = 0
sesgo_estimador_sigma_0 <- mean(resultados_0$p_gorro) - 0.5
print(sesgo_estimador_sigma_0)
## [1] -0.0006492864
# sigma = 1
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_1 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())

for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 1)
  
  # Almacenar los resultados
  resultados_1 <- rbind(resultados_1, data.frame(n = n, sigma = 1, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_1, aes(x = n, y = p_gorro)) +
  geom_line(color = "green") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 1)",
       x = "n",
       y = expression(hat(p)[n, 1])) +
  theme_minimal()

# Sesgo para sigma = 1
sesgo_estimador_sigma_1 <- mean(resultados_1$p_gorro) - 0.5
print(sesgo_estimador_sigma_1)
## [1] -0.03590337
# sigma = 2
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_2 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())

for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 2)
  
  # Almacenar los resultados
  resultados_2 <- rbind(resultados_2, data.frame(n = n, sigma = 2, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_2, aes(x = n, y = p_gorro)) +
  geom_line(color = "red") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 2)",
       x = "n",
       y = expression(hat(p)[n, 1])) +
  theme_minimal()

# Sesgo para sigma = 2
sesgo_estimador_sigma_2 <- mean(resultados_2$p_gorro) - 0.5
print(sesgo_estimador_sigma_2)
## [1] -0.02068251
# sigma = 4
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_4 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())

for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 4)
  
  # Almacenar los resultados
  resultados_4 <- rbind(resultados_4, data.frame(n = n, sigma = 4, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_4, aes(x = n, y = p_gorro)) +
  geom_line(color = "purple") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 4)",
       x = "n",
       y = expression(hat(p)[n, 1])) +
  theme_minimal()

# Sesgo para sigma = 4
sesgo_estimador_sigma_4 <- mean(resultados_4$p_gorro) - 0.5
print(sesgo_estimador_sigma_4)
## [1] -0.03427481

Respuesta

    • Como el sesgo es la diferencia entre el valor esperado del estimador y el verdadero valor del parámetro, al calcular el sesgo para los diversos valores de \(\sigma\), obtenemos que para \(\sigma=0\) y \(\sigma=1\) el sesgo es cercano a 0. El ruido añadido es muy pequeño o inexistente, por lo que el estimador está estimando correctamente el parámetro \(\hat{p}_{n,\sigma}\).
    • Para \(\sigma=2\), el sesgo es negativo, lo que sugiere que el estimador está subestimando el valor de \(\hat{p}_{n,\sigma}\).
    • Para \(\sigma=4\), el sesgo es positivo, lo que significa que el estimador está sobreestimando el valor de \(\hat{p}_{n,\sigma}\).
    • Por lo tanto los estimadores con menor sesgo tienden a ser más precisos, ya que en promedio sus valores están más cerca del verdadero valor de \(\hat{p}_{n,\sigma}\).

Pregunta 2: Intervalo de Confidencia

En las preguntas 2, 3 y 4 deberá trabajar con el dataset diabetes_prediction_dataset.

datos <- read.csv("diabetes_prediction_dataset.csv", header = TRUE)
head(datos)

Una forma de estimar la distribución de la media de una población es a través de la realización de la sampling distribution de una distribución cualquiera. El valor obtenido en cada sampling distribution nos entrega un estimador de la media que posee una determinada distribución de la población. Sabiendo esto, calcule la sampling distribution de las columnas bmi y blood_glucose_level, obteniendo el intervalo de confianza de \(95\%\) para cada una de las medias obtenidas desde la distribución. Para realizar este experimento considere los siguientes puntos:

Hints:

  • Para comparar comparar las 2 estimaciones puede ser útil normalizar bmi y blood_glucose_level a una misma escala.
  • Para realizar la sampling distribution podría serle útil el comando sample().
  • Puede ser útil generar un dataframe para graficar con ggplot2.

Respuesta:

# Definimos tamaños de muestreo
tamano_muestra = 100
n_muestras = 5000

# Inicializamos estructuras necesarias
list_mean_bmi = vector()
list_interval_bmi = vector()
list_typeCI_bmi = vector()
list_prop_bmi = vector()
sucess_bmi = 0

list_mean_bgl = vector()
list_interval_bgl = vector()
list_typeCI_bgl = vector()
list_prop_bgl = vector()
sucess_bgl = 0

#función para calcular desviación estándar
sd.p <- function(x){sd(x)*sqrt((length(x)-1)/length(x))}

# Obtenemos la media y desviación estándar de cada columna
mu_bmi <- mean(datos$bmi)
sd_bmi <- sd.p(datos$bmi)
mu_bmi
## [1] 27.32077
sd_bmi
## [1] 6.63675
mu_glucose <- mean(datos$blood_glucose_level)
sd_glucose <- sd.p(datos$blood_glucose_level)
mu_glucose
## [1] 138.0581
sd_glucose
## [1] 40.70793
alpha <- 0.05

# Sampling distribution, calculo del intervalo de confianza y proporción.

for (i in 1:n_muestras) {
  ### Para bmi ###
  samples_bmi <- sample(datos$bmi, tamano_muestra, replace = TRUE)
  sample_mu_bmi <- mean(samples_bmi)
  sample_sd_bmi <- sd(samples_bmi)
  
  # Calculo de intervalo de confianza
  se_bmi <- sample_sd_bmi/sqrt(tamano_muestra)
  error_bmi <- qnorm(1-alpha/2) * se_bmi
  left_bmi <- sample_mu_bmi - error_bmi
  right_bmi <- sample_mu_bmi + error_bmi
  
  # Se guardan los resultados 
  list_mean_bmi[i] <- sample_mu_bmi
  list_interval_bmi <- c(list_interval_bmi, left_bmi, right_bmi)
  list_prop_bmi[i] <- sucess_bmi/i
  
  # Se verifica si la media de la población se encuentra dentro del intérvalo de confianza
  if (left_bmi <= mu_bmi && right_bmi >= mu_bmi) {
    sucess_bmi <- sucess_bmi + 1
  }
  
  
  ### Para blood_glucose_level ###
  samples_glucose <- sample(datos$blood_glucose_level, tamano_muestra, replace = TRUE)
  sample_mu_glucose <- mean(samples_glucose)
  sample_sd_glucose <- sd(samples_glucose)
  
  # Calculo de intervalo de confianza
  se_glucose <- sample_sd_glucose / sqrt(tamano_muestra)
  error_glucose <- qnorm(1 - alpha / 2) * se_glucose
  left_glucose <- sample_mu_glucose - error_glucose
  right_glucose <- sample_mu_glucose + error_glucose
  
  # Se guardan los resultados
  list_mean_bgl[i] <- sample_mu_glucose
  list_interval_bgl <- c(list_interval_bgl, left_glucose, right_glucose)
  list_prop_bgl[i] <- sucess_bgl / i
  
  # Se verifica si la media de la población se encuentra dentro del intérvalo de confianza
  if (left_glucose <= mu_glucose && right_glucose >= mu_glucose) {
    sucess_bgl <- sucess_bgl + 1
  }
  
}

# Generamos dataframes para plotear mas facilmente los datos.

# Convertimos los vectores de intervalos en matrices de dos columnas
list_interval_bmi_matrix <- matrix(list_interval_bmi, ncol = 2, byrow = TRUE)
list_interval_bgl_matrix <- matrix(list_interval_bgl, ncol = 2, byrow = TRUE)

# Generamos dataframes para graficar más fácilmente los datos
df_bmi <- data.frame(media = list_mean_bmi, ci_low = list_interval_bmi_matrix[, 1], ci_high = list_interval_bmi_matrix[, 2], prop_ci = list_prop_bmi)
df_bgl <- data.frame(media = list_mean_bgl, ci_low = list_interval_bgl_matrix[, 1], ci_high = list_interval_bgl_matrix[, 2], prop_ci = list_prop_bgl)

# Verificar si el intervalo de confianza contiene la media de la población
df_bmi$contiene_media <- (df_bmi$ci_low <= mu_bmi & df_bmi$ci_high >= mu_bmi)
df_bgl$contiene_media <- (df_bgl$ci_low <= mu_glucose & df_bgl$ci_high >= mu_glucose)
# Plot de Intervalos de confianza
# Graficar una muestra aleatoria de 100 medias e intervalos de confianza para BMI
set.seed(123)  # Para reproducibilidad
sample_indices_bmi <- sample(1:n_muestras, 100)  # Seleccionar 100 índices aleatorios
sample_bmi <- df_bmi[sample_indices_bmi, ]

# Crear gráfico para la muestra de BMI
fp_bmi <- ggplot(data = sample_bmi, aes(x = factor(1:100), y = media, ymin = ci_low, ymax = ci_high)) +
  geom_pointrange() + 
  geom_hline(yintercept = mu_bmi, lty = 2) +  # Línea horizontal para la media poblacional
  coord_flip() +  # Girar coordenadas
  xlab("Muestra Aleatoria de Medias BMI") + ylab("Media (95% CI)") +
  theme_bw() + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  # Opcional: limpiar etiquetas del eje y
print(fp_bmi)

# Plot de Intervalos de confianza
# Graficar una muestra aleatoria de 100 medias e intervalos de confianza para Glucose
set.seed(123)  # Para reproducibilidad
sample_indices_bgl <- sample(1:n_muestras, 100)  # Seleccionar 100 índices aleatorios
sample_bgl <- df_bgl[sample_indices_bgl, ]

# Crear gráfico para la muestra de Glucose
fp_bgl <- ggplot(data = sample_bgl, aes(x = factor(1:100), y = media, ymin = ci_low, ymax = ci_high)) +
  geom_pointrange() + 
  geom_hline(yintercept = mu_glucose, lty = 2) +  # Línea horizontal para la media poblacional
  coord_flip() +  # Girar coordenadas
  xlab("Muestra Aleatoria de Medias Glucose") + ylab("Media (95% CI)") +
  theme_bw() + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  # Opcional: limpiar etiquetas del eje y
print(fp_bgl)

# Suponiendo que ya tenemos los datos de las muestras y las estadísticas
# Normalizamos las variables
normalizar <- function(x) {
  return((x - mean(x)) / sd(x))
}

# Normalizamos las medias obtenidas de las muestras
list_mean_bmi_normalized <- normalizar(list_mean_bmi)
list_mean_bgl_normalized <- normalizar(list_mean_bgl)

# Creamos un data frame para facilitar el plotting
df_comparacion <- data.frame(
  index = 1:n_muestras,
  bmi = list_mean_bmi_normalized,
  blood_glucose_level = list_mean_bgl_normalized
)

# Convertimos el data frame a formato largo sin reshape2
df_long <- data.frame(
  index = rep(df_comparacion$index, 2),
  Variable = rep(c("BMI", "Blood Glucose Level"), each = n_muestras),
  Valor = c(df_comparacion$bmi, df_comparacion$blood_glucose_level)
)

# Graficamos usando ggplot2
library(ggplot2)

ggplot(df_long, aes(x = index, y = Valor, color = Variable)) +
  geom_line() +
  labs(title = "Comparación de Estimaciones Normalizadas de BMI y Glucosa en Sangre",
       x = "Índice de Muestra",
       y = "Estimación Normalizada") +
  theme_minimal() +
  theme(legend.title = element_blank())

# Plot de proporción de CI
p_prop_bmi <- ggplot(df_bmi, aes(x = 1:n_muestras, y = prop_ci)) +
  geom_line(color = "blue") +
  ggtitle("Proporción de Intervalos de Confianza que Contienen la Media - BMI") +
  xlab("Iteraciones") + ylab("Proporción Acumulada") +
  ylim(0, 1)

p_prop_bgl <- ggplot(df_bgl, aes(x = 1:n_muestras, y = prop_ci)) +
  geom_line(color = "green") +
  ggtitle("Proporción de Intervalos de Confianza que Contienen la Media - Glucose") +
  xlab("Iteraciones") + ylab("Proporción Acumulada") +
  ylim(0, 1)

# Mostrar las gráficas
print(p_prop_bmi)

print(p_prop_bgl)

Respuesta Es posible observar que la estimación para blood_glucose_level tiene mejores resultados ya que a simple vista se observan más valores cercanos a la media real, en cambio para bmi se observan valores mas dispersos y no tan cercanos a la media real por lo que se espera eque este estimador tenga un mayor error.

Pregunta 3: Estimación de Máxima Verosimilitud

El objetivo de esta parte será obtener los parámetros de las distribuciones de la media y la mediana de bmi. Para responder la pregunta realice los siguientes puntos:

Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo la librería MASS).

Respuesta

# Obtenemos métricas de muestreo con repetición
set.seed(99)
means <- c()
medis <- c()
n_muestras <- 10000 
tamano_muestra <- 100 

for (i in 1:n_muestras) {
  muestra <- sample(datos$bmi, tamano_muestra, replace = TRUE)
  means <- c(means, mean(muestra))
  medis <- c(medis, median(muestra))
}

# dataframe con las medias y medianas de las muestras
resultados <- data.frame(media = means, mediana = medis)
# Histograma media
ggplot(resultados, aes(x = media)) + 
  geom_histogram(binwidth = 0.1, fill = "blue", color = "white") + 
  ggtitle("Histograma de Medias Simuladas") + 
  theme_minimal()

# Histograma mediana
ggplot(resultados, aes(x = mediana)) + 
  geom_histogram(binwidth = 0.1, fill = "green", color = "white") + 
  ggtitle("Histograma de Medianas Simuladas") + 
  theme_minimal()

# Media
# función log likelihood
ll_plot <- function(mu, sigma) {
  likelihood <- -sum(dnorm(means, mean = mu, sd = sigma, log = TRUE))
  return(likelihood)
}

ll_plot <- Vectorize(ll_plot)

# Definir el espacio donde se va a evaluar ll_plot
mu_m <- seq(26, 29, length.out = 10)    # Rango de valores para la media
sigma_m <- seq(0.01, 1, length.out = 10)   # Rango de valores para la desviación estándar

likelihood_m <- outer(X = mu_m, Y = sigma_m, ll_plot)

filled.contour(x = mu_m, y = sigma_m, z = likelihood_m, xlab = expression(mu), 
               ylab = expression(sigma), color = topo.colors, 
               main = "Log Likelihood para la Media")

# Mediana
ll_plot_mediana <- function(mu, sigma) {
  likelihood <- -sum(dnorm(medis, mean = mu, sd = sigma, log = TRUE))
  return(likelihood)
}
ll_plot_mediana <- Vectorize(ll_plot_mediana)
mu_d <- seq(27, 28, length.out = 10)
sigma_d <- seq(0.01, 1, length.out = 10)

likelihood_d <- outer(X = mu_d, Y = sigma_d, ll_plot_mediana)

filled.contour(x = mu_d, y = sigma_d, z = likelihood_d, xlab = expression(mu), 
               ylab = expression(sigma), color = topo.colors, 
               main = "Log Likelihood para la Mediana")

likelihood_mean <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  
  # Definimos la likelihood como la suma logarítmica de la función de densidad
  log_likelihood <- -sum(dnorm(means, mean = mu, sd = sigma, log = TRUE))
  
  return(log_likelihood)
}

# Agrrgue el rango donde operará nlminb
result_mean <- nlminb(objective = likelihood_mean, 
                      start = c(mu = 27, sigma = 1),  # Valores iniciales
                      lower = c(mu = 26, sigma = 0),  # Límite inferior
                      upper = c(mu = 28, sigma = 2))   # Límite superior
# Mostrar los resultados
cat("Estimación por Máxima Verosimilitud para la media:\n")
## Estimación por Máxima Verosimilitud para la media:
cat("Mu óptimo:", result_mean$par[1], "\n")
## Mu óptimo: 27.32434
cat("Sigma óptimo:", result_mean$par[2], "\n")
## Sigma óptimo: 0.6633704
likelihood_med <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  
  # Definimos la likelihood como la suma logarítmica de la función de densidad
  log_likelihood <- -sum(dnorm(medis, mean = mu, sd = sigma, log = TRUE))
  
  return(log_likelihood)
}

# Agrrgue el rango donde operará nlminb
result_med <- nlminb(objective = likelihood_med, 
                      start = c(mu = 27, sigma = 1),  # Valores iniciales
                      lower = c(mu = 26, sigma = 0),  # Límite inferior
                      upper = c(mu = 28, sigma = 2))   # Límite superior
cat("Estimación por Máxima Verosimilitud para la mediana:\n")
## Estimación por Máxima Verosimilitud para la mediana:
cat("Mu óptimo:", result_med$par[1], "\n")
## Mu óptimo: 27.30507
cat("Sigma óptimo:", result_med$par[2], "\n")
## Sigma óptimo: 0.1025047
# Muestra de medias
mu_media <- 27.32434
sigma_media <- 0.663372
n_muestra <- 10000
muestras_media <- rnorm(n_muestra, mean = mu_media, sd = sigma_media)
df_muestras_media <- data.frame(muestras_media)
# Histogrmas
hist_media_simulada <- ggplot(df_muestras_media, aes(x = muestras_media)) +
  geom_histogram(binwidth = 0.1, fill = "blue", color = "white") +
  ggtitle("Histograma de Medias Simuladas") +
  theme_minimal()
hist_media_simulada

# Muestra de medianas
mu_mediana <- 27.30507
sigma_mediana <- 0.1025046
n_muestra <- 10000
muestras_mediana <- rnorm(n_muestra, mean = mu_mediana, sd = sigma_mediana)
df_muestras_mediana <- data.frame(muestras_mediana)
# Histogrmas
hist_mediana_simulada <- ggplot(df_muestras_mediana, aes(x = muestras_mediana)) +
  geom_histogram(binwidth = 0.1, fill = "green", color = "white") +
  ggtitle("Histograma de Medianas Simuladas") +
  theme_minimal()
hist_mediana_simulada

Al analizar las distribuciones de medias y medianas obtenidas por máxima verosimilitud, se observa lo siguiente:

  • Medias: La estimación de máxima verosimilitud captura adecuadamente la variabilidad observada en los datos empíricos. Esto se refleja al comparar los histograma de las medias empíricas y el de las medias simuladas, lo que indica que la distribución normal estimada es un buen ajuste para los datos.

  • Medianas: En el caso de las medianas, no se observa el mismo comportamiento. El histograma de las medianas empíricas muestra una fuerte concentración en un solo valor, lo que indica una baja variabilidad en los datos originales. Esto difiere con las medianas simuladas, que presentan mayor variabilidad, propia de una distribución más normal.

La diferencia en la variabilidad entre las medianas empíricas y las simuladas sugiere que la distribución de los datos empíricos tiene características que no son bien capturadas por la estimación normal, posiblemente debido a una concentración excesiva de valores en un punto específico.

Pregunta 4: Bootstrap I

En esta pregunta se error e inetrvalo de confianza para la varianza de la columna HbA1c_level.

Las actividades por realizar son las siguientes:

Respuesta:

# Bootstrap
B <- 5000
largo <- length(datos$HbA1c_level)
output <- c()

for (it in 1:B) {
  sample_data <- sample(datos$HbA1c_level, largo, replace = TRUE)
  output[it] <- var(sample_data)
}

# Valor real de la varianza
real_variance <- var(datos$HbA1c_level)
real_variance
## [1] 1.146339
# Gráfico de puntos
plot(output, main = "Varianzas obtenidas con Bootstrap",
     ylab = "Varianza", xlab = "Muestra de Bootstrap")
abline(h = real_variance, col = "red", lty = 2) # Línea roja para la varianza real
legend("topright", legend = "Varianza real", col = "red", lty = 2)

# Histogrma
hist(output, main = "Histograma de varianzas de Bootstrap",
     xlab = "Varianza", breaks = 30, col = "lightblue")
abline(v = real_variance, col = "red", lty = 2) # Línea roja para la varianza real
legend("topright", legend = "Varianza real", col = "red", lty = 2)

# Obtenemos error e intervalos de confianze
se_vars <- sd(output)
z_a2 <- qnorm(1 - alpha/2)
alpha <- 0.05

# límite inferios
l.CI <- quantile(output, alpha/2)
# límite superior
u.CI <- quantile(output, 1-alpha/2)

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
## [1] "El intervalo de confidencia de 95% de las varianzas es: (1.136,1.157)"
sprintf('El SE de la varianzas es: (%.3f)', se_vars)
## [1] "El SE de la varianzas es: (0.005)"

Respuesta:

Para el gráfico de puntos para la varianza obtenida con boostrap se puede ver que los puntos están mayormente concentrados en la línea punteada roja, la cual corresponde a la varianza real. Además se nota que hay puntos outliers presentes pero en menor cantidad.

Para el histrograma de la varianza se cómo claramante las varianzas obtenidas con Bootstrap tienen una distribución cercana a la normal centrada en el valor de la varianza real.

Se puede inferir que al aplicar el método bootstrap para error estándar de la varianza, estos resultados tendrán una distribución cercana a la normal centrada en el valor de la varianza real, más aún en el histograma ya se puede apreciar cuál será un valor aproximado de el 95% de intervalo de confianza de la estimación.

Pregunta 5: Bootstrap II

El siguiente código genera una regresión lineal de bmi con respecto a age, es decir, una función lineal de age, \(y(age) = b + m\cdot age\), que pretende determinar el valor de bmi.

linearMod <- lm(bmi ~ age, data=datos)
print(linearMod)
## 
## Call:
## lm(formula = bmi ~ age, data = datos)
## 
## Coefficients:
## (Intercept)          age  
##    23.15536      0.09945
m <- linearMod$coefficients["age"]
b <- linearMod$coefficients["(Intercept)"]

ggplot() +
  geom_point(aes(x=datos$age, y=datos$bmi)) +
  geom_line(aes(x=datos$age, y=datos$age*m+b), color="red") +
  ggtitle("Regresión lineal") +
  ylab("bmi") + 
  xlab("Age")

Repita el proceso de la pregunta 4 para estimar el error e intervalos de confianza de los parámetros de la regresión (\(m\) y \(b\)). ¿Qué indican los resultados de Bootstrap?¿Un valor bajo en el error significa que la regresión es buena?

# Bootstrap
B <- 500
largo <- length(datos$bmi)
output_m <- numeric(B)
output_b <- numeric(B)

for (it in 1:B) {
  # Tomar una muestra con reemplazo de los datos
  sample_indices <- sample(1:largo, largo, replace = TRUE)
  sample_data <- datos[sample_indices, ]
  
  # Ajustar la regresión lineal con la muestra
  bootstrap_model <- lm(bmi ~ age, data = sample_data)
  
  # Guardar los coeficientes en las salidas correspondientes
  output_m[it] <- bootstrap_model$coefficients["age"]
  output_b[it] <- bootstrap_model$coefficients["(Intercept)"]
}

# Mostrar algunos resultados del bootstrap
head(output_m)
## [1] 0.09902192 0.09975384 0.10031168 0.10003177 0.10001440 0.09812545
head(output_b)
## [1] 23.16575 23.11355 23.10619 23.13401 23.11621 23.19340
# Obtenemos error e intervalos de confianze
se_m <- sd(output_m)
z_a2 <- qnorm(1 - alpha / 2)
alpha <- 0.05

# límite inferios
l.CI <- mean(output_m) - z_a2 * se_m
# límite superior
u.CI <- mean(output_m) + z_a2 * se_m

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
## [1] "El intervalo de confidencia de 95% de las varianzas es: (0.098,0.101)"
sprintf('El SE de la varianzas es: (%.3f)', se_m)
## [1] "El SE de la varianzas es: (0.001)"
# Obtenemos error e intervalos de confianze
se_b <- sd(output_m)
z_a2 <- qnorm(1 - alpha / 2)
alpha <- 0.05

# límite inferios
l.CI <- mean(output_b) - z_a2 * se_b
# límite superior
u.CI <- mean(output_b) + z_a2 * se_b

sprintf('El intervalo de confidencia de 95%% de los sesgos es: (%.3f,%.3f)', l.CI , u.CI)
## [1] "El intervalo de confidencia de 95% de los sesgos es: (23.156,23.159)"
sprintf('El SE del sesgo es: (%.3f)', se_b)
## [1] "El SE del sesgo es: (0.001)"

Respuesta

  • Error estándar bajo: Indica que los coeficientes son consistentes a través de las muestras, lo que sugiere un modelo de regresión estable.

  • Un error bajo no garantiza que el modelo sea bueno en términos predictivos, pero lo hace fiable de que observada relaciones en los datos actuales.

 


A work by CC6104